% randomly select N solution then local search
clear
%%
d = 2;
x_bound = repmat([0.25,10],[d,1]);
solutionId = fullfact(ones(1,d).*6);
OptSet.sol = exp((4.*solutionId-11)./20.*pi);
OptSet.fval = 0;
precision_init = 0.04; % 0.2912
OptSet.radius = 0.03;
MaximalBudget = 200000;
precision = [0.03, 0.0096, 0.003, 0.00096];

iepsilon = 1;
N = 40000;
radius = 0.2;
%% generate
X = rand([N,d]).*repmat(x_bound(:,2)'-x_bound(:,1)',[N,1])+repmat(x_bound(:,1)',[N,1]);
fitness = ObjValue(X);
%% ranking the solutions
[fitness, Ranking] = sort(fitness,'descend');
X = X(Ranking,:); % sort the samples from bad to good
SampleSet = [fitness,X];
%% clearing
OptimaSet = ones(N,1);
sampleid = 1;
while 1
    Distance = (sqrt(sum((repmat(X(sampleid,:),[N,1])-X).^2,2)))';
    neighbor = find(Distance<=radius);
    [~,best_neighbor] = min(fitness(neighbor));
    OptimaSet(neighbor) = 0; % delete all neighbor
    OptimaSet(neighbor(best_neighbor)) = 1; % save the best neighbor to the optima set
    next = find(OptimaSet((sampleid+1):end),1);
    if isempty(next)
        break
    else
        sampleid = sampleid + next;
    end
end
Optima_Stage1 = SampleSet(OptimaSet==1,:);
%% local search
step = sqrt(radius^2/d)/2;
step_threshold = precision(iepsilon)/4;
[Optima_Stage2, path, Budget_Stage2] = LocalSearch_BoxConstraint(@ObjValue,x_bound',Optima_Stage1,step,step_threshold);
sum(Budget_Stage2)
%% figure
if d==2
    % STAGE ONE
    figure()
    hold on
    scatter(X(:,1),X(:,2),2,fitness,'filled')
    plot(Optima_Stage1(:,2),Optima_Stage1(:,3),'x','MarkerSize',6,'LineWidth',2,'MarkerEdgeColor' ,[0.86,0.34,0.07])
    hold off
    % STAGE TWO
    figure()
    plot(Optima_Stage2(:,2),Optima_Stage2(:,3),'o')
end
